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Abstract 


A new Iterative scheme for solving boundary value problems Is 
presented. It consists of the Introduction of an artificial time depen- 
dence Into a modified version of the system of equations. Then explicit 
forward Integrations In time are followed by explicit Integrations 
backwards In time. The method converges under much more general condi- 
tions than schemes based In forward time Integrations (false transient 
schemes). In particular It can attain a steady state solution of an 
elliptical system of equations even If the solution Is unstable, In 
which case other Iterative schemes fail to converge. The simplicity of 
Its use makes It attractive for solving large systems of nonlinear 


equations . 



1. Introduction 


Elliptic equations are often solved numerically by Iterative 
methods. It is general'*}’ possible to associate an iterative scheme 
with some time marching scheme corresponding to a parabolic equation, 
l.e., the solution of the elliptic equation Is obtained as the steady 
state solution of a corresponding parabolic equation [1]. 

If the steady state solution of a parabolic system of equations 
Is unstable , developing Instabilities will prevent such solution from 
being reached starting from any set of :*nltlal conditions. In a similar 
way Instabilities will prevent "time marching" Iterative schemes from 
converging to an unstable solution. 

We have developed a new Iterative scheme with which It Is possible 
to find the steady state solution of a parabolic system of equations 
even If It Is unstable. The new scheme Is very simple to use and of 
rather general application. It can be used to solve boundary value 
problems, and, more generally, large systems of linear or nonlinear equa- 
tions not easily solved by other Iterative techniques. 

For simplicity the sche:ue is introduced first for the solution 
of a single complex equatlou (section ?) , and acceleration procedures 
are discussed In section 3. In section 4 the method Is extended to 
systems of equations. In section 5 we present examples of application 
to the solution of systems of equations, and section 6 contains a des- 
cription of three cases In which the method was used to obtain physically 
unstable steady state solutions. 
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2. Iterative solution of a single complex equation 


Consider the transcendental (complex In general) equation 


f(u) ■ 0. 


( 2 . 1 ) 


We Introduce an artificial time dependence and try to find the steady 
state solution of 


du 

dt 


f (u). 


( 2 . 2 ) 


If u ■ u Is a solution of (2.1) we write u u + v, and assuming v 
o o 

sufficiently small » expand (2.2) In a Taylor series 


^ - YV + 0(v2) . 


(2.3) 


Here y 



Is a complex number If f(u) Is complex. From 


u 

o 

(2.3) It Is clear that for any time marching scheme , there may be values 
of Y for which |v| will be amplified, so that Iterative schemes based on 
(2.2) will not converge. On the other hand the scheme that we will now 
present will converge to the solution for any value of y. 

Let us replace (2.2) by a modified equation 


du 

dt 


f*(u) , 


(2.4) 


where the asterisx represents the complex conjugate. The scheme consists 
of an explicit Euler "time" integration, followed by an explicit Euler 
Integration going backwards In "time". In both steps of the Iteration 
the modified equation (2.4) Is used: 

u ■ u'^ + At f (u'*) (2.5a) 

r 

u'^^ - 0 - At f*(a). 


(2.5b) 



We use the superindex v to denote Iteration number, not time, since we 
always return to the same "time level". 

When we linearize with respect to the departure from the solution, 
as In (2.3), we obtain 

^ ■ V + At Y v'^ (2.6a) 

vfl * * 

V ^ ^ - At Y ^ (2.6b) 

and eliminating 

« (1 - At^ yy )v'* . (2.7) 

This result Indicates that if the stability criterion • 

At^ yy* <2 (2.8) 


Is satisfied the time dependent component v will be damped out. This also 
shows that the method Is of first order convergence, l.e., the error Is 
reduced by an approximately constant factor at each Iteration so that It 

V . * 0 V V ' 

decreases as v 'u (1 - yy At^) . 

3f 

In the case of a double root, l.e., when v ■ — ■ 0 at u • u , 

’ ' 3u o 

. 3 v2 ^ 

the method can still be applied. Then f(u) •» (r^) + 0(v^), 

u«u ^ 

and when (2.5) Is applied we obtain 

V - (1 Y" nn V V )v (2.9) 


where n 
and v'^ 





at u « u . 
o 


Therefore the convergence rate Is much slower. 


The "time step" At Is a computational parameter, which may vary 
with V. For a linear equation It Is obvious that the optimum value Is 
At^P^ ■ . For a nonlinear equation y can be determined from the 
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equation Itself or from the solution after the first Iteration (see sec- 


tion 3), but for rapidly varying functions It Is safe to take Lt < dt^^^ 

and go "curve crawling" [2]. In practice, however, 

Lt can be determined as the largest value that produces stable results. 

Figure 1 shows graphically how the method works for the case of 

a real y. Values of for different values of u (assuming v constant) 

oc 

are plotted. Since the ~ Is proportional to “-u^, the derivatives 
are steeper for values of u far from the solution u^. If At ■ ■ -|^ 

the method converges In one Iteration (for linear equations). It should 
be noted that If, as In Figure 1, Re(y) > 0, a simple time marching scheme 


Is unstable. 


o* 


3. Acceleration of the scheme applied to a single equation 
3.1 Extrapolation 

When the error v la sufficiently small, y may be considered approxi- 
mately constant, and from the results of each Iteration we can extrapolate 
and obtain second order convergence. From (2.6) we get 

V * V 

u - a • -Y At V 


u'’ - 


v+1 


YY 


dtV 


(3.1) 


and therefore the best estimate of v 


V 


Is 




(u^-a)(u^-g)* 


. nd 


* 

Y At 



v+1 
u 

- a 


(3.2) 


(3.3) 


Then for the next Iteration optimum estimates of u and Lt are obtained 
as follows: 


vfl 
u . 
opt 


V _ (u^-Q) (u^-g)* 

. V v+1.* 
(u -U ) 


u + 0 (v^ ) 
o 


(3.4) 

✓ 


opt |y*V''1 


(3.5) 


3.2 Real equation 

If f(u) Is real, It Is necessary to apply only one of the forward 
and backward parts of scheme (2.5). In other words, (2.5) reduces to 

- o'* + At’f(u'') (3.6) 


where At Is determined from a single complete Iteration (2.5): 




tv '■ 



Fig. 2 indicates how v.he method works If this option is used. 
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4. Solution of a system of equations 

Consider the system of equations 

F(u) - 0 


(4.1) 


where F is a nonlinear, complex in general, matrix operator, and u is 
the vector of the dependent variables. If is a solution ve expand 
in a Taylor series with respect to the departure from the solution, and 
neglect terms higher than linear: 

F(u) « F(u ) + Jv - Jv . (4.2) 

o 

3F 

Here. J ■ I is the Jacobian such that ■ 3F /3u.. 

oU Ij 1. j 

U«U •' 

o 

Note that a simple extension of scheme (2.5) to higher dimension 


vi “ u'^ + At F (u'^) 
u'^ ■ fl - At F (a) 


(4.3) 


will not converge in general to the solution, because J can have complex 
eigenvalues even if F is real. 

Let J ■ S + A, where S ■ (J + J^)/2, A - (J - J^)/2 are symmetric 
and antisymmetric matrices respectively and the superindex T denotes 
transpose. Then the extension of scheme (2.6) to a system of equations is: 

* 

^ • v'^ + At (S + A ) v'^ (4.4a) 


« 

V ■ V 


V 


+ At (S 


* 

A ) 


(4.4b) 


.v+1 


- b - At 


. * * 

(S + A ) 


a* 

V 


(4.4c) 


Here, as before, the asterisk denotes complex conjugate. When b and v 
are eliminated in (4.4), we obtain 




'i '•;« 



n. 

- II - at2(S* + A*)(S - A)]v'' (4.5) 

where I It the Identity matrix. 

Therefore a sufficient condition for convergence la 

t 

S*k - A*S . (4.6) 

It Is easily verified that If J Is real, condition (4.6) Is equivalent 

T T 

to J being normal, l.e., JJ ■ J J. 

If condition (4.6) Is fulfilled, then from (4.5) and S ^ - S , 

St a 

A ■ -A we obtain 

v'^^ - [I - 4t2(S*^S + A*^A)]v''. (4.7) 

At at 

The matrix S S + A A Is the sum of positive definite matrices and thi*re- 
fore it Is positive definite, l.e.. It has real positive eigenvalues u^. 


If the linear stability criterion 

4t2p < 1 (4.8) 

n 

Is satisfied for all n, the method wllx converge. 

For a nonlinear problem, scheme (4.4) becomes ^ 

a - u'' + 4t(F*(u'') + F*(u'')) (4.9a) 

u u'' + At(F*(u'^) - F^(u'')) (4.9b) 

u'^^ ■ a - At(Fg(u) + (4.9c) 


where F_ and F are the nonlinear operators corresponding to S and A res- 
S A 

pectlvely. It should be noted that, computationally, the most expensive 
part of the scheme, which Is the evaluation of F(u), is executed only 
twice per Iteration, in steps (4.9a) and (4.9c). 
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It It possible to design an extrapolation procedure as in (3.2) 

y 

to determine v and obtain second order convergence. It would require, 
however, an explicit evaluation of the jacobian J and of its inverse, 
and therefore it would destroy the simplicity of the method. 


i 



12 . 


5. Examples of applications 

In the following examples we have not tried to optimize the procedure, 
but rather to teat the ability of f.he method to solve nonlinear equations. 

It should be noted that a good initial guess is nof necessary. 

5.1 Single transcendental equations 

a) Real equations 

i) We found roots of u ■ tan u b> determining the uly 
state solutions of ■ (u - tan u)/u. We divided by u in order to avoid 

0 L 

the trivial solution u ■ 0. The conservative value At ■ .2 At . was 

opt 

used. Once two succeslve approximations differed by less than 1% the 
*extrapr1 ation procedure was used. lO-AO Iterations were needed to get 
a va'.u.^ c:ose to the solution, depending on the initial guess, and then 
1>3 extrapolations were enough to converge with an error smaller than 
10 When the option described in section 3.2 was used, about one 
fourth the number of iterations were ne'cessary. 

ii) sin u ■ .5. Using the option described in section 3.2 
only 2-4 Iterations and 1-2 extrapolations are necessary for convergence. 

The value At ■ 1 was determined from the equation Itself. 

ill) sin u ■ 1. This is a double root case, and convergence 
is very slow. It takes about 100 iterations to reduce the error to 
IX. 

b) Complex equations 

We solved sin u • 1.17520 using the standard scheme (2.3). 

Starting from Re(u) ■ 2, lm(u) ■ 0.5, At ■ 1, 10 iterations and 2 
extrapolations yielded the solution Re(u^) ■ 1.570796, Im(u^) - 0.583628. 

The same result was obtained after 25 standard iterations. 
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5.2 Solution of ordinary differential equations 

a) Steady state solution of ■ (u+2)~ + (u+1) , with periodic 

o C oX 

boundary conditions: u(0) u(l). The only steady state solution Is 

■ -1, but It Is unstable. We found Iteratively the steady state 
solution by applying (4.9) In the following way: 

u «= u'^ + At { (u'* + 2)u^ + (o'* + 1)} 

u = u'^ + At {-(u'^ + 2)u^ + (u'^ + 1)} 

u'^^ « fl - At{ (u'^ + 2)u^ + (u + 1)}. 

The x-derlvatlves on the right hand side were replaced by centered 

o 

finite differences. We started from the Initial guess u = sin 2inc 

and used a time step At ® .1. About 400 Iterations were necessary to 

obtain the solution u ■ -1 with an error smaller than IZ. 

o 

b) A nonlinear boundary value problem 

This method can be used to solve nonlinear boundary value pro- 
blems In a rather Inefficient way, but with minimum human effort. 
Consider the example given by Acton [2], page 171: 

y" " y^(l “ 0.2 sin 2x) * 0 

y(0) = 1, y(l) = 0 . 

Since there are no first derivatives the problem Is solved In 
the following way, using finite differences to express x-derlvatlves: 

y “ y'* + 4t{y^ - (y'^)2(l - 0.2 sin 2x)} 

y'^1 . y - At{y^ - (y)2(l - 0.2 sin 2x)} . 

About 300 Iterations were necessary to converge within 1%, and after 
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800 Iterations the error was smaller than 10 
c) A "stiff" equation 

The equation - 101 + lOOy ■ 0 has solutions ■ e* and 

y^ “ e^^^. Only y^^ satisfies the boundary conditions y(0) ■ 1, 
y(l) “ e, but the equation cannot be solved by shooting or by relaxa- 
tion because the "hidden solution" y^ is amplified. 

We tried to solve this equation by writing 


y - y'* + At{y^ - lOly^ + lOOy'*] 

(5.1a) 

y - y'* + At{y^ + lOly^ + lOOy''} 

(5.1b) 

^ - 4t{y^ - lOly^^ + lOOy} 

(5.1c) 


and replacing the x-derivatlves by centered second order finite differences 
in the interval (0,1). After about 5000 iterations the solution remained 
bounded but it failed to converge. A closer analysis shows that the scheme sat- 
isfies condition (4.5) everywhere except at the interior points next to 
the boundaries, where the finite differences corresponding to the first 
and second derivatives do not commute. This is easily solved by 
eliminating y, y in (5.1) and then applying finite differences on 


v+1 

y 


y'"- LtHy^ 

^ •'xxxx 


- lOOOly^ + lOOOOy'^} 


at the points next to the boundaries. With 20 grid intervals and 

-4 -4 

At ■ 2 X 10 the method converged with an error smaller than 10 

in about 2000 iterations. 


5.3 A highly nonlinear system of equations 

Scarf [3] developed a constructive method to find fixed points 
of a mapping that transforms a space into itself. He tested the method 
on a highly nonlinear system of equations 
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-i . .. 


£ ( 

i-i 


k»l \ 
i a 

i k-1 Ak \ 


-wj 


1 - 1,...,5 


(5.2) 


5 

subject to the normalization condition ir^ ■ 1. Here a^^, and 
are positive data matrices, and ir^ Is a vector of unknowns. According 
to Scarf [3] this system of equations cannot be solved by gradient 
methods. 

We solved the problem by multiplying (5.2) by and applying 

only 4.9a and 4.9c under the as amptlon F^(ti) * 0* The solution was 
normalized after each Iteration. The method converged In 300 Iterations 
within an error of 10 However when we tried a simple forward scheme 
(4.9a), It also converged and In only 60 Iterations. This Illustrates 
the point that whenever a simple forward scheme converges. It should 
be used, because It Is more efficient than the forward-backward scheme 
(cf. section 3.2). 

For comparison we note that Scarf's [3] constructive algorithm 
converged In 168 Iterations, but the averaged solution given by Scarf 
has errors varying between 2 and 20Z. 


6. Physically unstable steady state solutlvons 

As discussed In the Introduction, this scheme was originally 
devised as a method of finding steady state solutions of partial differen- 
tial equations which are unstable, so that small perturbations are ampli- 
fied by regular Iterative techniques. 

The method has been tested and found to yield an unstable steady 
state solution In three different cases. 

6.1 Rotating annulus: example of a spectral (Fourier) model 

Lorenz [4] developed a simplified spectral model of a rotating 
annulus heated In Its rim and cooled In the center. It consists of a 
system of only 14 ordinary differential equations for the real and 
Imaginary Fourier components. For certain values of the external para- 
meters, the thermal Rossby number and the Taylor number, the steady 
state solution becomes unstable and Rossby waves develop. 

In a spectral (Fourier) model the linearized Jacobian J (see 
eq. 4.2) Is complex, but It has only diagonal elements. Therefore the 
straightforward scheme (4.3) can be used. 

Following Lorenz [4] we made a regular forecast for an unstable 
situation, until Rossby waves were well developed. Then, using the same 
time step as Lorenz, we applied (4.3) and, after 300 Iterations, the 
model converged to the unstable steady state solution found analytically 
by Lorenz. 

6.2 An unstable steady state boundary layer 

Inez Fung [5] has used scheme (4.9) to obtain the unstable 
steady state axlsymmetrlc nonlinear boundary layer of a geophysical vor- 
tex, where the Coriolis parameter Is different from zero. Two cases 
were considered: a) a V-R vortex, and b) Carrier et al. [6] model 


hurricane vortex. Convergence occurred In about 1000 iterations in both 
cases, and the results for case b) compared well with those for the 
downdraft region In Carrier's model. 

6.3 A steady state solution with shear instability: example of a finite 

differences model 

The flow defined by the equation 

1^ - + + C sin y 

has a steady state solution 

C ^ 

5 » — sin y 

o V ^ 

which is unstable due to the presence of inflection points in the pro- 
file of the mean velocity field. Here C is the vorticlty, the stream 
function Ip Is defined through the equation 

- c , 

V Is a diffusion coefficient and C Is a constant. 

We used centered finite differences and periodic boundary condi- 
tions to solve the problem. Starting from the analytical steady state 
solution we made a regular forecast, using Lorenz' [7] N-cycle scheme 
with N » 8. Immediately strong Instabilities developed, until the mean 

flow was completely distorted. Then we started again from the analytical 

C 

solution but adding a strong perturbation of the form sin 2x sin 2y, 
and applied scheme (A. 9). After about 800 iterations the analytical 


solution was recovered. 
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7. Summary and conclusions 

We have presented a new numerical scheme to solve nonlinear 
boundary value problems. It consists of the Introduction of an arti- 
ficial time dependence Into a modified versior. of the equations. Then 
explicit forward Integrations In "time" are followed by explicit Inte- 
grations backwards In "time". For a real system of equations. If only 
the forward part of the step Is executed, the method reduces to one of 
the class of "false transient" methods that has been widely used to 
solve elliptical systems of equations. (See discussion by Varga, [1] ^ 

pp. 270-278, and Malllnson and de Vahl Davis, [8] ) The forward- 
backward method Is less efficient than the forward methods, since It 
requires twice as many computations per Iteration and the convergence 
rate Is usually slower. On the other hand, the forward-backward scheme 
converges under much more general conditions than the forward schemes. 

In particular. It converges to the steady state solution of an ellip- 
tical system of equations even If the solution Is unstable. In which 
case the forward schemes fall. 

The method has been tested for several cases, and the numerical 
examples Indicate that when the size of the time step is chosen close 
to the maximum value compatible with numerical stability, convergence 
occurs on the order of serveral hundred Iterations. 

The number of Iterations necessary for convergence Is highly 

dependent on the size of the time step, and on the difference between 

*X *T 

largest and smallest eigenvalues of the matrix S S + A A (section 4). 
Therefore a thoughtful variation of the time step size with the itera- 
tion number, as done for exa.aple in the Okamura-Rivas scheme (Grant and 
Rivas, [9] may increase significantly the efficiency of the method. 

It may also be possible to normalize the system of equations in order 

P.TTPX'r'' ^ ' 'I 
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to decrease the range of the eigenvalues This is equivalent to 

choosing different time scales for different equations. 

The numerical method presented here, though strictly iterative, 
shares all the simplicity of Initial value problems solved with expli- 
cit time differencing schemes. It Is hoped that It will prove useful 
for solving systems of large numbers of nonlinear equations not easily 
solved with other Iterative techniques. 
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Figure Captions 


Fig. 1 
Fig. 2 


: Graphic Interpretation of the Iterative scheme (2.5) for a real 

equation f(u) ■■ 0. y Is assumed constant and the dashed lines 
Indicate the time derivative for different values of u. 

: Same as Fig. 1, except for the scheme (3.6). 
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